/*
 * Solves the Panfilov model using an explicit numerical scheme.
 * Based on code orginally provided by Xing Cai, Simula Research Laboratory
 *
 * Modified and  restructured by Scott B. Baden, UCSD
 */

#include <assert.h>
#include <stdlib.h>
#include <iostream>
#include <iomanip>
#include <string.h>
#include <math.h>
#include <sys/time.h>
#include <mpi.h>

using namespace std;
#define BUFSIZE 2048


// Utilities
//

// Timer
// Make successive calls and take a difference to get the elapsed time.
static const double kMicro = 1.0e-6;
double getTime()
{
	struct timeval TV;
	struct timezone TZ;

	const int RC = gettimeofday(&TV, &TZ);
	if(RC == -1) {
		cerr << "ERROR: Bad call to gettimeofday" << endl;
		return(-1);
	}

	return( ((double)TV.tv_sec) + kMicro * ((double)TV.tv_usec) );

}  // end getTime()

// Allocate a 2D array
double **alloc2D(int m,int n){
	double **E;
	int nx=n+1, ny=m+1;
	E = (double**)malloc(sizeof(double*)*ny + sizeof(double)*nx*ny);
	assert(E);
	int j;
	for(j=0;j<ny;j++) E[j] = (double*)(E+ny) + j*nx;
	return(E);
}

// Reports statistics about the computation
// These values should not vary (except to within roundoff)
// when we use different numbers of  processes to solve the problem
double stats(double **E, int m, int n, double *_mx){
	double mx = -1;
	double l2norm = 0;
	int i, j;
	for (j=1; j<=m+1; j++)
		for (i=1; i<=n+1; i++) {
			l2norm += E[j][i]*E[j][i];
			if (E[j][i] > mx)
				mx = E[j][i];
		}
	*_mx = mx;
	l2norm /= (double) ((m+1)*(n+1));
	l2norm = sqrt(l2norm);
	return l2norm;
}

// External functions
extern "C" {
	void splot(double **E, double T, int niter, int m, int n);
}
void cmdLine(int argc, char *argv[], double& T, int& n, int& px, int& py, int& do_stats, int& plot_freq, int& no_comm);

void append(char *str, double x)
{
	char buf[10];
	sprintf(buf, "%5.2f",x);
	strcat(str,buf);
}

// Main program
int main (int argc, char** argv)
{
	/*
	 *  Solution arrays
	 *   E is the "Excitation" variable, a voltage
	 *   R is the "Recovery" variable
	 *   E_prev is the Excitation variable for the previous timestep,
	 *      and is used in time integration
	 */
	double **E, **R, **E_prev;

	// Various constants - these definitions shouldn't change
	const double a=0.1, b=0.1, kk=8.0, M1= 0.07, M2=0.3, epsilon=0.01, d=5e-5;

	double T=1000.0;
	int m=200,n=200;
	int do_stats = 0;
	int plot_freq = 0;
	int px = 1, py = 1;
	int no_comm = 0;

	cmdLine( argc, argv, T, n,px, py, do_stats,  plot_freq, no_comm);
	m = n;

	/*START CODE ADDED*/

	//+2 to take the weird indexing into account
	int xsize = m+2;
	int ysize = n+2;
	//the dimensions of each process array
	int xinterval = (m+1)/px;
	int yinterval = (n+1)/py;
	//    printf("proc num: %d", procnum);

	//cout << "py is: " << py << endl;

	cout << "py is: " << py << endl;
	cout << "px is: " << px << endl;
	cout << "x interval is " << xinterval << " m is " << m <<endl;
	cout << "y interval is " << yinterval << " n is  " << n << endl;

	// Allocate contiguous memory for solution arrays
	// The computational box is defined on [1:m+1,1:n+1]
	// We pad the arrays in order to facilitate differencing on the
	// boundaries of the computation box
	E = alloc2D(xinterval+3,yinterval+3);
	E_prev = alloc2D(xinterval+3,yinterval+3);
	R = alloc2D(xinterval+3,yinterval+3);
	//should it be x,ysize +1?
	double **allE = alloc2D(xsize-1,ysize-1);

	int i,j;

	//No MPI functions called before this
	MPI_Init(&argc, &argv);


	int procnum;
	int procrank;
	//get the number of processes
	MPI_Comm_size(MPI_COMM_WORLD, &procnum);
	//get the rank of the process
	MPI_Comm_rank(MPI_COMM_WORLD, &procrank);


	//initialize the process array
	//horizontal decomposition
	int ylowbound;
	int yhighbound;
	int xlowbound;
	int xhighbound;
	printf("proc num: %d", procnum);
	cout << "py is: " << py << endl;
	cout << "px is: " << px << endl;
	cout << "x interval is " << xinterval << " m is " << m <<endl;
	cout << "y interval is " << yinterval << " n is  " << n << endl;

	ylowbound = 1;
	xlowbound = 1;
	if (procnum == 1) {
		no_comm = 1;
		xhighbound = n+1;
		yhighbound = m+1;
	} else if (px == 1) {
		yhighbound = yinterval+1;
		xhighbound = n+2;
		//vertical decomposition
	} else if (py==1) {
		xhighbound = xinterval+1;
		yhighbound = m+2;
	}

	// Initialization:
	//if vertical decomposition
	for (j=1; j<=yhighbound; j++)
	{
		for (i=0; i<=xhighbound; i++)
		{
			E_prev[j][i] = R[j][i] = 0;
		}
	}

	printf("xlowbound: %4d, proc: %2d\n", xlowbound, procrank);
	printf("ylowbound: %4d, proc: %2d\n", ylowbound, procrank);
	printf("xhighbound: %3d, proc: %2d\n", xhighbound, procrank);
	printf("yhighbound: %3d, proc: %2d\n", yhighbound, procrank);

	if (procnum == 1)
	{
		// serial
		printf("serial decomp\n");
		for (j=1; j<=m+1; j++)
		{
			for (i=n/2+2; i<=n+1; i++)
			{
				E_prev[i][j] = 1.0;
			}
		}
		for (j=m/2+2; j<=m+1; j++)
		{
			for (i=1; i<=n+1;i++)
			{
				R[j][i] = 1.0;
			}
		}

	}
	//horizontal decomposition
	else if(px ==1)
	{
		
//               printf("horizontal decomp\n");
//		cout << "process #"<<procrank<<endl;
                for (j=ylowbound; j<=yhighbound-1; j++)
                {
		cout << "process #"<<procrank<< " for j as " << j<<endl;

                        for (i=(xhighbound+1)/2; i<=xhighbound-1; i++)
                        {
//				cout << "process #"<<procrank << " for j as "<<j<<" for i as "<<i<<endl;

                                E_prev[i][j] = 1.0;
                        }
                }
		cout << "process #"<<procrank<<" SECOND TIME"<<endl;
                //check if over halfway y-point of display indices
                if(procrank >= procnum/2)
                        for(i=xlowbound; i <= xhighbound-1;i++)
                                for (j=ylowbound; j<=yhighbound-1; j++)
                                        R[i][j] = 1.0;

	//vertical decomposition
	} else if(py == 1){
		printf("vertical decomp\n");
		//check if over halfway x-point on display indices
		if(procrank >= procnum/2)
			for (j=xlowbound; j<=xhighbound; j++)
				//should the low-bound here be +1 or +2?
				for (i=(yhighbound-2)/2+2; i<=yhighbound; i++)
					E_prev[j][i] = 1.0;

		for(j = (xhighbound-1)/2+2; j <= xhighbound;j++)
			for (i=ylowbound; i<=yhighbound; i++)
				R[j][i] = 1.0;

	}

	//DEBUG PRINT STATEMENTS.............
/*
	cout << "E_prev for process: "<< procrank <<" ylowbound " << ylowbound-1 << " yhighbound " << yhighbound << " xlowbound-1 " << xlowbound-1 << "xhighbound " << xhighbound <<endl;
	for (j=ylowbound-1; j<=yhighbound; j++)
	{
		printf("proc: %2d, j: %3d --> ", procrank, j);
		for (i=xlowbound-1; i<=xhighbound; i++){
			//cout << " j:" << j << " i:" << i << "  "<<endl;
			printf("%1.0f",E_prev[i][j]);
		}
		cout << endl;

	}

	cout << "E for process: "<< procrank<< endl;
	for (j=ylowbound-1; j<=yhighbound; j++) {
		printf("proc: %2d, j: %3d --> ", procrank, j);
		for (i=xlowbound-1; i<=xhighbound; i++){
			//cout << " j:" << j << " i:" << i << "  "<< E[j][i];
			printf("%1.0f",E[i][j]);
			//cout << E[i][j];0 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0

		}
		cout << endl;
	}
	cout << "R for process: "<< procrank<<endl;
	for (j=ylowbound-1; j<=yhighbound; j++)
	{

		printf("proc: %2d, j: %3d --> ", procrank, j);
		for (i=xlowbound-1; i<=xhighbound; i++){
			//cout << " j:" << j << " i:" << i << "  "<< R[j][i];
			printf("%1.0f",R[i][j]);
		}
		cout <<  endl;
	}

	cout << endl;*/
	/*END CODE ADDED*/

	double dx = 1.0/n;


	// For time integration, these values shouldn't change
	double rp= kk*(b+1)*(b+1)/4;
	double dte=(dx*dx)/(d*4+((dx*dx))*(rp+kk));
	double dtr=1/(epsilon+((M1/M2)*rp));
	double dt = (dte<dtr) ? 0.95*dte : 0.95*dtr;
	double alpha = d*dt/(dx*dx);

	cout << "dx is: " << dx << endl;

	cout << "dt= " << dt << ", T = " << T << endl;
	cout << "processor geometry: " << px << " x " << py << endl;


	if (no_comm)
		cout << " *** Communication DISABLED" << endl;

	cout << endl;

	// Start the timer
	double t0 = MPI_Wtime();


	/*CODE ADDED*/


	//number of borders of each process
	int numborders;
	//if top or bottom, only 1 border, otherwise 2 borders
	if(procrank == 0 || procrank == procnum-1)
		numborders = 1;
	else
		numborders = 2;


	//MPI state variables
	MPI_Status* status = (MPI_Status*)malloc(numborders*sizeof(MPI_Status));
	MPI_Request* request = (MPI_Request*)malloc(numborders*sizeof(MPI_Request));

	/*END CODE ADDED*/

	// Simulated time is different from the integer timestep number
	// Simulated time
	double t = 0.0;
	// Integer timestep number
	int niter=0;



	char str[BUFSIZE];
	sprintf(str,"proc: %d, BEFORE MESSAGE PASSING\n", procrank);
	for (j=ylowbound-1; j<=yhighbound; j++)
	{
		for (i=xlowbound-1; i<=xhighbound; i++){
			append(str,E_prev[i][j]);
		}
		strcat(str, "\n");

	}
	printf(str);


	int count = 0;


	while (t<T) {

		t += dt;
		niter++;

		/*CODE ADDED*/
		//no idea why we need a tag?
		int tag = 0;



		//use message passing here to communicate the border regions to each process
		if(!no_comm) {
//			cout << "reaches here 0 " << procrank << endl;

			//the upper border
			if(procrank > 0){
				printf("proc: %d, do upper border12345\n", procrank);
				char str[BUFSIZE];
				double recvCells[xinterval];
				double sendCells[xinterval];


//				cout << "reaches here 1 " << procrank << endl;
				// get the ghost borders and ptu into the cells array
				MPI_Irecv(&recvCells[0], xinterval, MPI_DOUBLE, procrank-1,
						MPI_ANY_TAG, MPI_COMM_WORLD, request );


//				cout << "reaches here 2 "<< procrank  << endl;
				
                                for(i=1;i<=xinterval;i++){
                                        sendCells[i-1] = E_prev[i][1];
                                }

				sprintf(str,"proc: %d, Sending to %d:\n", procrank, procrank-1);
				for(i=0;i<xinterval;i++){
					append(str, sendCells[i]);
					//printf("%2.2f ", E_prev[xhighboun1][i]);
				}
				strcat(str, "\n");
				printf("sending:\n");
				printf(str);
//				cout << "reaches here 3 "<< procrank  << endl;


				// fill up the array of cells to send
				//send border region to subsequent process (top most row)
				MPI_Send(&sendCells[0], xinterval, MPI_DOUBLE, procrank-1, 0, MPI_COMM_WORLD);
//				cout << "reaches here 4 "<< procrank  << endl;



				//wait for recieving...
				MPI_Wait(request,status);
//				cout << "reaches here 5 "<< procrank  << endl;

				count++;
				// fill in the ghost borders that we received
				for (i=1;i<=xinterval;i++) {
					//E_prev[i][0] = recvCells[i-1];
					E_prev[i][0] = recvCells[i-1];

				}
//				cout << "reaches here 6 " << procrank << endl;

				sprintf(str,"proc: %d, Received from %d:\n", procrank, procrank -1);
				for(i=0;i<xinterval;i++){
					append(str,recvCells[i]);
					//printf("%2.2f ", E_prev[xhighbound][i]);
				}
				strcat(str, "\n");
				printf("recieved: \n");
				printf(str);
//				cout << "reaches here 1 " << procrank << endl;

			}

			//the lower border
			if(procrank < procnum-1)
			{
				printf("proc: %d, do lower border12345\n", procrank);
				
				double recvCells[xinterval];
				double sendCells[xinterval];
				char str[BUFSIZE];
				//printf("%d: has south border\n", procrank);
				//accept the border region from preceeding process
				MPI_Irecv(&recvCells[0], xinterval, MPI_DOUBLE, procrank+1,
						MPI_ANY_TAG, MPI_COMM_WORLD, request );



				// fill up the array of cells to send
				for(i=1;i<=xinterval;i++){
					//printf("i: %d, xhb: %d, xint: %d\n", i, xhighbound, xinterval);
					//sendCells[i-1] = E_prev[i][xhighbound-1];
					sendCells[i-1] = E_prev[i][yhighbound-1];
				}


				sprintf(str,"proc: %d, Sending to %d:\n", procrank, procrank + 1);
				for(i=0;i<xinterval;i++){
					append(str,sendCells[i]);
				}
				strcat(str, "\n");
				printf("sending: \n");
				printf(str);

				//send border region to subsequent process (bottom most row)
				MPI_Send(&sendCells[0], xinterval, MPI_DOUBLE, procrank+1, 0, MPI_COMM_WORLD);

				//wait for recieving...
				MPI_Wait(request,status);
				count++;
				// fill in the ghost borders that we received
				for (i=0;i<xinterval;i++)
				{
					//E_prev[i][yhighbound] = recvCells[i-1];
					E_prev[i+1][yhighbound] = recvCells[i];
				}

				sprintf(str,"proc: %d, Receive from %d:\n", procrank, procrank + 1);
				for(i=0;i<xinterval;i++){
					append(str,recvCells[i]);
				}
				strcat(str, "\n");
				printf("recieved: \n");
				printf(str);
			}
		}
		/*END OF CODE ADDED*/

		cout << "MADE IT TO HERE " << endl;
		/*
		 * Copy data from boundary of the computational box
		 * to the padding region, set up for differencing
		 * on the boundary of the computational box
		 *
		 */
		if (px == 1)
		{
			//printf("fill in ghost cells for horizontal decomp\n");

			// if there is no north border, then go ahead and replace the
			//	upper ghost cells
			if (procrank == 0)
			{
				for (j=1; j<= xinterval; j++) {
					E_prev[j][0] = E_prev[j][2];
				}
			}

			// replace the left and right ghost cells, respectively
			for (i=1; i<=yinterval; i++)
			{
				E_prev[0][i] = E_prev[2][i];
			}

			for (i=1; i<=yinterval; i++)
			{
				E_prev[xhighbound][i] = E_prev[xhighbound-2][i];
			}


			// if there is no south border, then go ahead and replace the
			//	lower ghost cells
			if(procrank == procnum-1)
			{
				for (j=1; j<=xinterval; j++) {
					E_prev[j][yhighbound] = E_prev[j][yhighbound-2];
				}
			}


/*			printf("proc: %d, ghost cells===========================\n", procrank);
			for(i=0;i<=xhighbound;i++){
			for(j=0;j<=yhighbound;j++){
				cout<<E_prev[i][j];
			}
				cout<<endl;
			}
*/

			// Solve for the excitation, the PDE
			for (j=1; j<xhighbound; j++){
#pragma ivdep
				for (i=1; i<yhighbound; i++) {

					E[j][i] = E_prev[j][i]+alpha*(E_prev[j][i+1]+E_prev[j][i-1]-4*E_prev[j][i]+E_prev[j+1][i]+E_prev[j-1][i]);
					//printf("proc %d, E[j][i]: %5.2f\n", procrank, E[j][i]);
/*					printf("procrank: %d, "
							"\tE_prev[%d][%d]: %5.2f\n "
							"\talpha: %21.10f\n"
							"\tE_prev[%d][%d]: %5.2f\n"
							"\tE_prev[%d][%d]: %5.2f\n"
							"\tE_prev[%d][%d]: %5.2f\n"
							"\tE_prev[%d][%d]: %5.2f\n", procrank,j, i,E_prev[j][i], alpha, j, i+1,E_prev[j][i+1],
							               j, i-1,E_prev[j][i-1],j+1,i, E_prev[j+1][i],j-1,i,E_prev[j-1][i]);*/
				}
			}


		}

		char str[BUFSIZE];
		sprintf(str,"proc: %d, AFTER MESSAGE PASSING %d\n", procrank, count);
		for (j=ylowbound-1; j<=yhighbound; j++)
		{
			for (i=xlowbound-1; i<=xhighbound; i++){
				//cout << " j:" << j << " i:" << i << "  "<< E_prev[j][i];
				append(str,E_prev[i][j]);
			}
			strcat(str, "\n");
		}
		printf(str);

		/*
		 * Solve the ODE, advancing excitation and recovery to the
		 *     next timtestep
		 */
		for (j=1; j<=xinterval; j++){
			double *RR = &R[j][1];
			double *EE = &E[j][1];


			for (i=1; i<=yinterval; i++, EE++,RR++)
			{
				EE[0] += -dt*(kk*EE[0]*(EE[0]-a)*(EE[0]-1)+EE[0]*RR[0]);
			}
		}

		for (j=1; j<=xinterval; j++){
			double *RR = &R[j][1];
			double *EE = &E[j][1];
			for (i=1; i<=n+1; i++, EE++, RR++)
				RR[0] += dt*(epsilon+M1* RR[0]/( EE[0]+M2))*(-RR[0]-kk*EE[0]*(EE[0]-b-1));

		}

//		double **tmp = E; E = E_prev; E_prev = tmp;

		//get all processes finished w/ computing before proceeding
		if(!no_comm) {
			MPI_Barrier(MPI_COMM_WORLD);
		}
		//   k = (int)(t/100);
		//   int k = (int)(t/10);
		if (do_stats){
			//       if ((t-k*100)<dt)
			if ((niter %  1000) == 0 ){
				cout << procrank << ": iteration " << niter << ", t = " << t << endl;
			}
		}
		if (plot_freq){
			//cout << " IN THE PLOT FREQUENCY PART" << endl;

			int k = (int)(t/plot_freq);
			if(no_comm){
				if ((t-k*plot_freq)<dt)
					splot(E,t,niter,m+1,n+1);
			}else{
			  //horizontal decomposition
			  if(px==1){	 
			     cout << "allE's dimensions are: x --> " << xsize-1 << " y --> " << ysize-1 << endl;
			     if((t-k*plot_freq)<dt){
				MPI_Barrier(MPI_COMM_WORLD);
				double **tmp_rcv = alloc2D(xinterval,yinterval);
			
				for(i=1; i<procnum; i++){
					if(procrank == 0){
						for(i = 0; i < yinterval; i++)
							for(j = 0; j < xinterval; j++)
								allE[j][i] = E[j+1][i+1];

						MPI_Status status;
						MPI_Irecv(&tmp_rcv[0][0], (xinterval*yinterval), MPI_DOUBLE, MPI_ANY_SOURCE, 0,MPI_COMM_WORLD, request);
						MPI_Wait(request, &status);
						int source = status.MPI_SOURCE;
						int offset_start = (source)*yinterval;
						int offset_end = (source+1)*yinterval;
						cout << "offset start " << offset_start << " offset end " << offset_end << endl;
						cout << "source: " << source<<" yinterval/procnum "<<(yinterval/procnum)<<endl;
						for(i = 0; i < yinterval; i++){
							for(j = 0; j < xinterval; j++){
								allE[j][offset_start+i] = tmp_rcv[j][i];
							}
						}
					}else{
						//horizontal decomposition
						if(px == 1){
							double sendCells[xinterval*yinterval];
							int counter = 0;
							cout << "xinterval " << xinterval << " y interval " << yinterval<< endl;
							for(i = 1; i<=yinterval; i++){
								for(j=1;j<=xinterval;j++){
									sendCells[counter] = E[j][i];
									cout << sendCells[counter]<< " ";
									counter++; 
								}
							      cout << endl;
							}
							MPI_Send(&sendCells[0], (yinterval*xinterval), MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
						}					
					}
				}
				//plot the data after having gathered all the data
				if(procrank == 0){
					if ((t-k*plot_freq)<dt){
						char str[BUFSIZE];
						sprintf(str,"proc: %d, allE array \n", procrank);
						for (i=0; i<ysize-1; i++)
						{
							for (j=0; j<xsize-1; j++){
							append(str,E[j][i]);
							}
							strcat(str, "\n");
						}
						printf(str);
						splot(allE,t,niter,(xsize),(ysize));
					}
				}
			    }
			}
		   }
		}//if plotfreq bracket
			//sleep(0.5);
			//       scanf("%d",&resp);
	
	double **tmp = E; E = E_prev; E_prev = tmp;

     }
	double t1 = MPI_Wtime();	
	double time = t1-t0;
	double final_time = 0.0;
	MPI_Reduce(&time, &final_time, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

	//	MPI_Barrier(MPI_COMM_WORLD);
	if(procrank == 0){
		
		cout <<  "End at time: " << T << ",  iteration " << niter << endl;
		double mx;

		double l2norm = stats(E_prev,xinterval,yinterval,&mx);
		printf("Max: %13.6e, L2norm: %13.6e\n",mx,l2norm);
		cout << "Wallclock time: " << final_time << " sec." << endl;
		if (plot_freq)
		{
			cout << "\n\nEnter any input to close the program and the plot..." << endl;
			int resp;
			scanf("%d",&resp);
		}
	}
	MPI_Finalize();
	/*END CODE ADDED*/
	free (E);
	free (E_prev);
	free (R);
}
